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We demonstrate the new class of variance reduction techniques for hadron propagator and nu- 
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varieties of observables compared to existing techniques such as low-mode averaging (LMA). We 
adopt this technique to hadron two-point functions and three-point functions, and compare with 
LMA and traditional source-shift method in the same ensembles. We observe AMA is much more 
cost effective in reducing statistical error for these observables. 
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1. Introduction 

In order to precisely evaluate non-perturbative quantities in lattice calculation, the reducing 
noise-to-signal ratio is one of the most important tasks especially for nucleon electric dipole mo- 
ment [jl|, the hadronic contribution to muon anomalous magnetic moment |J, nucleon form 
factors and structure functions [|5j|, X]' meson mass and mixing angle [^] and so on. We consider 
new strategies to effectively increase statistics without generating the new gauge configurations. 

Traditionally translational symmetry on the lattice has been used to increase the statistics of 
correlation function (correlator) of hadron interpolating operator(s). Since correlators for different 
source locations with the same relative separations between operators are exactly same in the infi- 
nite statistics limit for translational invariant action, the average over several source locations can 
be regarded as several times more statistical samples if correlation among correlators with differ- 
ent source is negligible. In this case, however, the additional computation of conjugate gradient 
(CG) at each source locations is needed. Low-mode-averaging (LMA) [0, |, ||] takes advantage 
of low-(eigen)mode for its dominant observables e.g. pseudoscalar correlator [f7j |J ^j, |RJ II, [120, 



to remove the cost of additional computation of CG. However for the other observables, nucleon 



propagator or resonance state and heavy mesons, [13, 14 15] reported that the statistical error 
reduction in LMA is less significant than pseudoscalar correlator. 

Here we examine the different variance reduction techniques from LMA for nucleon corre- 



lator and its three-point function. All-mode averaging (AM A) proposed in [16] reduces variance 
from a//-modes, without introducing bias, and thus should serve as a powerful technique to pre- 
cisely evaluate the observables including highly composite contribution coming from a//-modes. 
In this proceedings, we explain this idea and show the numerical results of hadron spectroscopy 
and isovector form factors of nucleon in realistic lattice setup. 

2. Covariant approximation averaging 

The observable G is calculated over gauge ensemble {U\ ,U2,--- , Un coo1 } to obtain the ensem- 
ble average; 

1 Nconl /i 

^) = j^Lm} + o(N; f), (2.i) 

where the second term denotes the uncertainty due to finite number of available gauge configura- 
tions. Under a transformation g G G, where G is a set of symmetry transformations on the lattice, 
the link variable is transformed as U — > U g . Ensemble average of an observable G should be 
covariant, 

(G*[U}) = (G[US}), (2.2) 

in the infinite statistics N con f — > °°. G 8 is the transformed observables. 

Here we introduce the approximation G^ appx ^ which fulfills the following condition: 
Appx-1: 

(AGAG^) ((AG) 2 )~{(AG^) 2 ), (2.3) 

^((AG) 2 ){(AG^)) 2 ) 

with AX = X - (X). r is the correlation between G and G^ ppx \ 

Appx-2: The computational cost of G^ avpx ^ is much smaller than original G. 
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Appx-3: Covariance of approximation : {0 a w[U g }) = (^p^s [£/]}. 
Using <^( a PP x ) we construct improved estimator; 

^.(imp) = ^.(rest) + ^ a PP x ) ; ^(rest) _ & _ ^(appx^ ^ a PP x ) = N^ 1 V ^( a PP x )«. (2.4) 

geG 

The first and the second conditions are to reduce statistical error of <^( rest ) and the computational 
cost respectively at fixed N con f, while the third one is to avoid the bias : it leads to (<^( im P)) = (&). 
When we perform Nq times measurements for <^( a PP x ) after transformation g G G, for instance, 
shifting source locations, the statistical error of improved estimator, A( imp ), will be reduced to 



A (imp) 2(1 -r)+N G \ (2.5) 

compared with original error A ignoring the correlation among ^( a PP x )§ with different g's. In the 

-1/2 

case of strong correlation, r ~ 1, in Appx-1, A( imp ) becomes nearly N G times smaller than A. 
If Cost(^ appx) ) is significantly cheaper than N G xCost(^) (Appx-2), total cost of 0^ for a 
fixed size of error is reduced by a factor ~ Cost (^ appx) ) /Cost (^) + N G l (the last factor comes 
from computation of 6 in <^( rest )). The above improved estimator defines covariant approximation 
averaging ( CAA ). 

LMA [§,|9|] is one of CAA; <£?( a PP x ) consists of low-mode, and g is a shift of source location. In 
LMA, <^( a PP x ) is correlator constructed from a few low-modes of Hermitian Dirac operator H{x,y) 
(or its even-odd preconditioned counterpart), where we only present formula for the point source 
case for simplicity, 

S^\x,y) = £ V'lW^M, 4^ = 7T I ^(S (low)g ), (2.6) 
k=\ ly ° geG 

with eigenmode Yk and eigenvalue Xu in Y*yH{x,y)Wk{y) = ^kYk( x )- Appx-3 is satisfied since 
Yk (x) respects the translational symmetry. Cost of ^( a PP x) includes the I/O for eigenmode, outer 



products of eigenvectors in Eq^ljS), and the cost to contract quark propagators for &(S^ ow ^), which 
is typically small cost (Appx-2). Appx-2,3 are satisfied for most of observables in LMA, however, 
Appx-1 strongly depends on G and size of Nx ~ 0(100). 

We proposed another example of CAA, which we call as all-mode averaging (AMA) in [|T6|]. 



Using the sloppy CG Q17[ ] combined with low-mode deflation (e.g flllQ) in which the stopping 
condition e of CG is made loose as £ama < 10~ 3 -10~ 4 (or fix the number of CG iterations to some 
small number l ), the approximation is given by 

S (all) (*,30 = I V'yfctoYi W +MB(x,y)), 4 AMA) = — Z ff(S^), (2.7) 
k=l ly ° geG 



'With the fixed stopping condition of CG in the approximation, S' a "), a bias will be introduced due to the finite 
precision (64 bits arithmetic in our case) breaking Appx-3 in theoretically very small provability. This bias could be 
avoided by fixing the iteration number to a constant as pointed out by M.Luscher and S. Hashimoto independently. In 
this proceedings, however, we checked that this bias is undetectable by comparing results of AMA and LMA, latter of 
which does not have such bias: for pion propagator they are identical within 1-cx (0.1%) error, as shown in Figure [lj 
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Table 1: Parameters in LMA/AMA. Ranges of CG iteration numbers in each ensembles are shown. 



m 


-Wconf 


N G 


Nx 


£ 


CG iter. 


£ama 


CG iter. (AM A) 


0.005 


380 


32 


400 


10" 8 


350-360 


3 x 10~ 3 


70-90 


0.01 


257 


32 


180 


io- 8 


600-630 


3 x 10~ 3 


90-130 



where f e denotes the polynomial function of H implicitly created in CG process to approximate 
inverse : / e (A) ~ l/X for Xn < A < A raax - AM A has advantage that S^ an \x,y) takes account of not 
only low-mode contribution but also (approximately) all-mode contribution which is controlled by 
the two parameters and £ama- AMA also fulfills the above three conditions (Appx-l-Appx-3) 
for a much wider class of observables than LMA. 

3. Numerical results 

We use the Nf = 2 + 1 domain-wall fermion (DWF) configurations generated by RBC/UKQCD 
collaboration in 24 3 x64 lattice at j3 = 2.13 Iwasaki gauge action [|l8|]. CG algorithm with four 
dimensional even-odd preconditioning was used to compute quark propagators at quark mass 
m = 0.01,0.005, and 5th dimension size is L s = 16. To calculate eigenmode of Hermitian even-odd 
preconditioned kernel of DWF operator in 10 _8 accuracy, we implement the implicitly restarted 
Lanczos algorithm with Chebychev polynomial acceleration Jl9|]. Note that in use of even-odd 
bases, one needs to choose the four dimensional shift vector of source point to meet Appx-3 (even 
steps in four directions are sufficient). 

In this proceedings, LMA/AMA estimator is obtained from ^ appx ^ with Nc = 32 different 
source locations separated by every 12 for spatial direction and 16 for temporal direction; (0,0,0,0), 
(12,0,0,0), (12,12,0,0),- • • ,(12,12,12,48) in lattice unit. (0,0,0,0) is the original source location for 
G. Stopping condition e for original observable and £ama for the sloppy CG in AMA, are 
defined as \\Hx — b\\/\\b\ \ < £,£ama with the even-site source vector b and the even-site solution 
vector x (see also Table [j]). Note that the number of CG iterations is for the case of deflated CG 
with Nx low-modes. To compare the performance, we set the Gaussian-type smearing source with 
parameters taken from [^J. In Ref.[^|] they measured three- and two-point functions for four source 
locations in temporal direction to extract the nucleon isovector form factor (and axial charge), and 
thus 4 x A^ C onf samples were accumulated. For m = 0.005, they have shown the results with quark 
sources set in two timeslices per one CG (double source method), by which effective samples are 
doubled. Furthermore their nucleon source was non-relativistic definition (2 spinor sources rather 
than 4). In our cost analysis below, we will correct these two factors for fair comparison. 

3.1 Two-point function 

In Figure [TJ, we compare nucleon (AO, pseudoscalar (P) and vector (V) meson correlator on 
three different time separations for original and LMA/AMA analysis. As mentioned before, LMA 
suppresses fluctuation from low-modes of Dirac matrix in Eq. ( $L(\ ) by averaging over No source 
locations of low-modes, and thus the error reduction of LMA is significant for larger time separa- 
tions. On the other hand, since AMA approximately averages contribution from all modes, AMA 
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Figure 1: The comparison between LMA/AMA and original analysis for nucleon (N), pseudoscalar (P) 
and vector (V) meson propagator at different time-slices t = 4,8, 12. The colored bar indicates the ratio of 
relative error between original and LMA/AMA. This is the case at m = 0.005. 



Table 2: The nucleon mass obtained by global fitting of N correlator with Gaussian smeared sink. We 
use GeV unit. Last two columns are the gain for AMA with and without deflation, which are inverse of 
reductions of costs to achieve the same size of errors. For Gain w / def , the cost for computing eigenvectors 
are ignored while Gain w / odef reflects its cost, see [[kJ for more detail breakup of cots. Costs of TY [||| are 
estimated in case of the relativistic source on a single time-slice per a CG for the fair comparison; 1864 
meas. (m — 0.005) and 1424 meas. in m = 0.01. 





ra/vCOrig) 


m N (LMA) 


m N (AMA) 


m N (TY) 


Gain w / def 


Gain w / odef 


0.33 


1.124(22) 


1.145(9) 


1.139(4) 


1.148(10) 


4 


15 


0.42 


1.221(17) 


1.219(11) 


1.233(4) 


1.217(9) 


4 


5 



is expected to suppress the fluctuation of observables in both short and long distances. The above 
expectation is clearly seen in Figure [I]; as a time-separation changes from t = 12 to t = 4, the error 
reduction ratio (green bar) is degraded in LMA significantly, while, for AMA, error reduction rate 

1 /2 

remains close to the ideal magnitude, I/Nq ~ 0.18, and always smaller than that of LMA, for 
every channels and distances. Error reduction for P channel is less different between LMA and 
AMA. This is likely because the single pion is mostly dominated by the low-modes, also discussed 
in ©111. 

In Figure ||, the effective mass of nucleon in AMA becomes flatter for both point and Gaussian 
smeared sinks compared to those of LMA. In Table || we see that the precision of nucleon mass in 
AMA is higher than previous study [||] while the computational cost is roughly less by 1/4 times. 
The detailed comparison between them including computational time of low-mode is discussed in 

Ph. 
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Figure 2: Effective mass plot of nucleon correlator in original, LMA/AMA analysis with point sink and 
Gaussian smeared sink. The colored bound shows the statistical error when globally fitting the propagator. 
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Figure 3: Isovector form factor F\ (q 2 ) and F2 (q 2 ) obtained in LMA/AMA and presented in TY et al.[||] at 
m = 0.005. The bars denote the relative error for each q 2 results. 



3.2 Nucleon isovector form factor 

To test AMA for more involved observables, we compare the three-point functions from 
AMA/LMA with parameters shown in Table [l] and the previous result in Nucleon isovec- 
tor form factors, which are particularly feasible test bed, are extracted from the ratio of three and 
two point functions: 



RJt\,t,to\p\,po) =K 



■'o,0) 



(h - 1 ,q)C N G (t - tp , 0)C^ (h - 1 , 0) 

c N L (h 



.1/2 



tQ,q)Cl(h-to,q) 



(3.1) 



Here C^ G (t,q) is the two-point function of nucleon with local (L) or Gaussian (G) sink at spatial 



momentum q, Cj (q, t) is the three-point functions with vector current , and K = ^/ 2{E^ + m^) / . 
See [|5|] for more details. Isovector form factors extracted from Eq.(|3.1[) are shown in Figure |[ 
Precision of F\^{(r) in AMA is better than the previous results [||] at all transfer momenta q 2 
examined. This demonstrates that AMA is effective to reduce errors for many lattice observables. 
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4. Summary 

In this proceedings, we have shown several results using the new class of (bias-free) error 
reduction techniques called as covariant approximation averaging (CAA). We proposed all-mode- 
averaging (AMA), in which contributions from all eigen modes are taken into account in the aver- 
aged approximation with the sloppy CG using deflation. AMA is applicable to a broad variety of 
observables including nucleon spectrum, three-point functions and other composite correlators. We 
compared the nucleon mass and isovector form factors for realistic Nf = 2+1 DWF configurations 
with lattice volume (2.7 fm 3 ) and light quark masses (m n ~ 0.3-0.4 GeV). Significant gain in com- 
putational cost, up to 15 times is observed when compared with traditional many-source method 
without the deflation (and 4 times compared with the case with deflation). AMA significantly re- 
duces error for cases in which LMA can not so much. Using this technique, calculations of the 
nucleon electric dipole moment and the hadronic contributions to the muon g-2 are now underway. 

Numerical calculations were performed using the RICC at RIKEN and the Ds cluster at FNAL. 
This work was supported by the Japanese Ministry of Education Grant-in-Aid, Nos. 22540301 
(TI), 23105714 (ES), 23105715 (TI) and U.S. DOE grants DE-AC02-98CH10886 (TI) and DE- 
FG02-92ER40716 (TB). We also thank BNL, the RIKEN BNL Research Center, and USQCD for 
providing resources necessary for completion of this work. 
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